npj Systems Biology and Applications

Anticancer drug resistance is challenging to overcome because it can arise through both intrinsic and acquired mechanisms, each driven by distinct cellular machinery. In particular, there is a sharp need for therapies that target hormone-insensitive prostate tumors due to the growing incidence of castration-resistant prostate cancer. Optimizing the pathways that regulate apoptosis in prostate cancer offers a promising strategy to induce apoptosis and inhibit tumor progression, since these mechanisms do not depend on hormonal signaling. Here, we identified strategies to enhance apoptosis in prostate cancer cells. We used several computational tools (including sensitivity analysis, particle swarm optimization, and ImageJ) to design an ordinary differential equation model of caspase-mediated prostate cancer apoptosis signaling. We apply the model to identify key modalities that increase the propensity toward apoptosis across three separate pro-apoptotic drugs (Tocopheryloxybutyrate, Narciclasine, and Celecoxib). Overall, we demonstrate that apoptosis dynamics can be accurately captured in response to each of the three drugs and identify which features of the model represent viable targets for overcoming intrinsic drug resistance.

AbstractsO_ST_ABSBackgroundC_ST_ABSQuantitative systems pharmacology (QSP) models provide mechanistic insight into drug response but are limited by labor-intensive, expert-driven workflows. We developed an AI-assisted QSP (AI-QSP) framework that integrates large language models (LLMs) with SBML-based modeling to enable automated reconstruction, extension, and calibration of mechanistic models. MethodsThe framework was applied to a published CAR-T QSP model. The model was reconstructed in SBML and extended via LLM-guided prompts to incorporate key resistance mechanisms: T-cell exhaustion, PD-1/PD-L1 checkpoint regulation, and tumor antigen escape. Model development followed an iterative expert-in-the-loop workflow. The resulting model (21 reactions, 9 species) was calibrated to synthetic benchmark data using 19-parameter optimization. Model credibility was assessed using ASME V&V 40 and ICH M15 principles, including global sensitivity and profile-likelihood analyses. ResultsThe calibrated model reproduced benchmark dynamics with high accuracy (mean log-RMSE = 0.132). Sensitivity analysis identified CAR-T killing and bystander cytotoxicity as dominant drivers of tumor response. Profile-likelihood analysis showed 71% of parameters were practically identifiable, with remaining parameters prioritised for future data-driven refinement. ConclusionsAI-assisted QSP modeling enables reproducible, scalable model reconstruction and evolution while maintaining mechanistic transparency and regulatory alignment. This framework provides a foundation for accelerating model-informed drug development in cell and gene therapies.

Experimental methods for characterizing single cells and cell populations have improved tremendously over the past decades. This progress has enabled the development of quantitative, mechanistic models for cellular processes based on either single cell or bulk data. However, coherent statistical frameworks for the model-based integration of different data types at the single-cell and population levels are still missing. In this work, we present a mathematical modeling approach for integrating single-cell time-lapse, single-cell snapshot, single-cell time-to-event and population-average data. Utilizing a formulation based on nonlinear mixed-effect modeling, we enable the description of multiple data types, with and without single-cell resolution, and we propose a tailored parameter estimation method. Furthermore, we propose a tailored parameter estimation scheme that facilitates the assessment of underlying process parameters. Our study demonstrates that the proposed approach can reliably integrate diverse data types, thereby improving parameter identifiability and prediction accuracy. Applying this framework of extrinsic apoptosis reveals that simultaneously considering multiple data types can be essential, particularly when experimental constraints limit data availability. The proposed approach is broadly applicable and may significantly advance our understanding of complex biological processes.

Circadian rhythms are self-sustained biological oscillations that coordinate diverse physiological processes in plants, including growth, metabolism, and environmental responses. These rhythms arise from an interconnected transcriptional translational feedback network that integrates multiple entrainment cues such as light and temperature. The plant circadian clock is organized around key regulatory loops involving CCA1, LHY, PRRs, TOC1, ELF4, LUX, and other transcriptional regulators, whose coordinated interactions ensure precise and robust oscillations. In this study, we developed an ordinary differential equation based mathematical model, building upon a previous framework to incorporate additional regulatory modules and transcriptional controls that better reflect experimentally observed behaviour. To elucidate the regulatory organization of this model, we performed a multi-layered computational analysis combining four complementary approaches: (i) period sensitivity analysis to quantify how parameter perturbations influence the systems timing, (ii) phase portrait analysis to visualize dynamic interactions among key components, (iii) knockout analysis to identify parameters essential for sustained rhythmicity, and (iv) network impact analysis using composite weighted network indices to evaluate hierarchical control across the network. Together, these analyses reveal that transcriptional repression, protein degradation, and light-regulated synthesis form the dominant control mechanisms within the circadian system. The results highlight a hierarchical and robust network structure centred on the CCA1/LHY and PRRs feedback loop, with redundant modules ensuring stability under perturbations. Thus, this model provides an improved, biologically consistent framework for dissecting the dynamic architecture of the plant circadian clock and guiding future experimental validation.

Alzheimers disease (AD) is a multifactorial neurodegenerative disorder characterized by coordinated dysregulation of multiple genes, requiring system-level approaches to elucidate underlying molecular mechanisms. While existing computational studies largely focus on differential expression analysis or machine-learning-based feature selection, they often overlook inter-gene relationships and network topology, limiting biological interpretability. In this study, we present a network-based framework for prioritizing candidate genes in Alzheimers disease by integrating gene co-expression network analysis with multiple centrality measures. Transcriptomic data comprising approximately 39,000 genes across 324 Alzheimers and control samples were preprocessed using log-transformation, variance filtering and Z-score normalization, followed by LASSO-based feature selection to retain phenotype-associated genes. A weighted gene co-expression network was constructed using Pearson correlation to capture coordinated transcriptional activity. Network topology was analyzed using degree, betweenness and eigenvector centrality to identify genes that are highly connected, act as information brokers or occupy influential positions within the network. A consensus ranking was derived by merging these centrality measures, enabling robust prioritization of candidate genes. The results highlight a subset of highly central genes, including several small nucleolar RNAs and regulatory transcripts implicated in RNA processing, synaptic function and neurodegenerative pathways. By jointly leveraging co-expression structure and complementary centrality metrics, the proposed framework provides an interpretable and reproducible strategy for identifying biologically meaningful gene candidates, offering potential utility for biomarker discovery and therapeutic target prioritization in Alzheimers disease.

Recent studies have established that the circadian clock influences onset, progression and therapeutic outcomes in a number of chronic conditions including cardio-metabolic diseases and cancer. For the latter, they also suggest that chronotherapy offers the potential to refine current treatments and improve the development of future anti-cancer medicines. Therefore, there is a need for tools to measure the functional state of the tumoural circadian clock in patients. We have previously developed a model-led machine-learning algorithm called TimeTeller which has the potential to provide such a tool. Here we demonstrate its potential for clinical relevance by a study of 1286 breast cancer patients in which we characterise the nature of the circadian clock disruption in their tumours and demonstrate a strong nonlinear association between 10-year survival and TimeTellers tumoural clock disfunction score {Theta}. We find that good tumour clock function is antagonistic to survival.

Cancer cells in hypoxic environments often proliferate less but exhibit enhanced migration relative to their normoxic counterparts. Recent in vitro and in silico studies have characterized the role of hypoxic memory - the ability of cancer cells to retain their hypoxic phenotype even when reoxygenated - in tumor invasion. However, the observations have been limited either to exposing cancer cells to hypoxia for a fixed duration or by assuming a fixed-time persistence of the hypoxic state upon reoxygenation independent of the duration of hypoxia exposure. Thus, time-dependent cell-state changes during hypoxia and their impact on hypoxic memory remains unclear. Here, we first analyze transcriptomic data from breast cancer samples to show that the genes upregulated at transcriptional level and hypomethylated at epigenetic level are enriched in cell invasion, indicating hypoxic memory-driven process of tumor invasion. Next, we used a computational model to investigate how the spatial-temporal dynamics of oxygen levels in a tumor drive time-dependent changes in hypoxic memory and influence tumor invasion dynamics. Our simulation results show that such dynamic hypoxic memory can drive enhanced tumor invasion over a fixed hypoxic memory by a) enriching hypoxic cell density at the tumor front, b) reducing sensitivity of hypoxic cell state to fluctuations in oxygen supply, and c) enhancing effective diffusion of hypoxic cells. Our results highlight the crucial role of dynamic hypoxic memory in shaping tumor invasion dynamics, underscoring the need to elucidate its underlying mechanisms in future studies.

Predicting cellular behaviors, a central task in systems biology and metabolic engineering, can be enhanced through integrative modeling of processes such as gene regulation and metabolism. Information flow from gene regulation (modeled via a gene regulatory network) to metabolism (modeled via a genome-scale metabolic model) is well-studied, but the reciprocal regulation of genes by metabolites is less explored. We introduce CausalFlux, a method that models bidirectional feedback between genes and metabolites, in order to predict steady-state reaction fluxes under wild-type (WT) or perturbed (e.g., gene knockout/KO) conditions. CausalFlux does so by iteratively performing causal surgery on a Bayesian gene regulatory network and constraint-based analysis of a coupled metabolic model. CausalFlux enabled us to assess the impact of two-way feedback in several testbed models and real-world biological systems by comparing its predictions to those of TRIMER, a state-of-the-art model of gene-to-metabolite one-way feedback. Incorporating bidirectional feedback, as in CausalFlux, improved the Spearman correlation between actual and predicted fluxes in 92% of the 39 distinct simulation conditions relative to TRIMER. For predicting growth/no-growth phenotype following single-gene KOs in E. coli, CausalFlux achieved a balanced accuracy of 0.79 in identifying essential genes, and TRIMER achieved 0.71 for the same task, again highlighting the importance of modeling two-way feedback. In ablation studies that further dissect the role of specific metabolite[->]gene feedback edges in E. coli, the F1 scores of gene essentiality predictions decreased by 7.5% and 13% upon ablation of feedback edges from any metabolite to the crp gene and the 10 metabolic feedback genes with the highest influence on the KO genes, respectively. Finally, we highlight the application of CausalFlux to predict the essentiality of several hundred genes under different media conditions. Overall, our findings show that CausalFlux can crucially utilize information on feedback metabolites to predict trends in reaction fluxes and qualitative (growth/no-growth) outcomes; thereby encouraging future systems modeling efforts to carefully incorporate not only gene-to-metabolite but also metabolite-to-gene interactions. AvailabilityCode pertaining to the CausalFlux method, and its benchmarking and application is publicly available at: https://github.com/BIRDSgroup/CausalFlux. Author summaryThe myriad processes within a living cell, such as gene regulation or metabolism, are tightly interconnected. Modeling these interconnected processes can offer a deeper mechanistic understanding of cellular behaviors, as well as guide efforts that engineer the metabolic output of a cell. In this work, we develop a novel integrated model of gene regulation and metabolism that incorporates bidirectional feedback between these two processes, via the concept of metabolite-induced causal surgery on a gene regulatory network and gene-induced constraints on the fluxes of metabolic reactions. Our model, which we call CausalFlux, represents an advance over most existing models that capture just the one-way gene-to-metabolism feedback (i.e., genes coding for enzymes that control metabolic reactions). Our CausalFlux methodology opens up an unique opportunity to quantify the impact of two-way feedback in gene-metabolite systems, via comparison of CausalFluxs predictions to those of TRIMER, a published model incorporating one-way feedback alone. For predicting reaction fluxes in testbed models and essential genes in E. coli, quantitative comparison of the performance of CausalFlux vs. TRIMER showed that accounting for two-way feedback leads to more accurate and biologically meaningful predictions. CausalFlux also enabled us to quantify the effect of two-way feedback by comparing prediction performance before and after ablation of certain feedback edges from metabolites to genes. Overall, our findings highlight the importance of modeling gene regulation and metabolism as two-way interconnected systems within a living cell, and encourage future works to incorporate gene{leftrightarrow}metabolite feedback into their analyses.

Constraint-based metabolic modeling is a powerful way to study the mechanistic basis of cellular states and disease, but effective use demands substantial computational expertise and careful coordination of multi-step analyses. We developed MechAInistic to lower this barrier enabling researchers to ask complex biological questions in natural language. MechAInistic is a multi-agent system harnessing large language models organized around an Architect-Reviewer pattern that that converts a natural-language question into an executable, model-grounded workflow and produces a structured report. It supports pathway comparison, perturbation analysis, drug-target exploration, and literature interpretation across healthy and disease paired states. We evaluated MechAInistics therapeutic hypothesis generation using two immune-cell use-cases. For rheumatoid arthritis/healthy Naive B models, it identified mitochondrial metabolic rewiring and nominated Devimistat/CPI-613 as an investigational OGDH-centered hypothesis. In CD4+ Th17 multiple sclerosis/healthy models, the workflow identified NADP-dependent isocitrate dehydrogenase as the optimal target and proposed Ivosidenib as an FDA-approved repurposing candidate. GRAPHICAL ABSTRACT O_FIG O_LINKSMALLFIG WIDTH=200 HEIGHT=83 SRC="FIGDIR/small/723319v1_ufig1.gif" ALT="Figure 1"> View larger version (19K): org.highwire.dtl.DTLVardef@1b5c1d1org.highwire.dtl.DTLVardef@1c798cforg.highwire.dtl.DTLVardef@10161d3org.highwire.dtl.DTLVardef@1bd7dce_HPS_FORMAT_FIGEXP M_FIG C_FIG

Metabolic networks are typically viewed as homeostatic systems that stabilize flux, energy charge, redox balance, and metabolite availability under perturbation. However, it remains unclear whether the same feedback architectures that support metabolic robustness can also generate learning-like, experience-dependent adaptation. Here, we develop a coarse-grained dynamical model of mammalian energy metabolism to test whether prior perturbation can improve future metabolic responses. The model represents core glucose, glutamine, fatty acid, and oxidative phosphorylation pathways as coupled ordinary differential equations with Michaelis-Menten-type fluxes, product-inhibition feedback, adaptive enzyme-capacity regulation, and explicit ATP costs for enzyme adjustment. Rather than aiming to reproduce quantitative fluxes for a specific cell type, the framework is designed to expose how metabolic feedback, regulatory cost, repeated perturbation, and environmental variability interact. We use this model to ask whether adaptive enzyme regulation enables improved recovery after repeated challenges, whether such effects depend on energetic control costs, and whether environmental variability broadens or constrains the set of reachable adaptive states. This approach provides a tractable way to investigate how homeostatic metabolic regulation may give rise to experience-dependent metabolic plasticity.

Over a thousand fungal genomes have been sequenced, yet manually curated genome-scale metabolic models (GEMs) are available for only a limited number of species. Moreover, these models have often been developed independently, leading to inconsistencies in namespaces, compartment definitions, and pathway representations that hinder comparative analysis, the systematic reuse of prior curation efforts, and the integration of consolidated metabolic knowledge. Here, we present the Consolidated Fungal Core Metabolism Model (CFCMM), constructed by integrating thirteen published fungal models spanning Ascomycota, Mucoromycota, and both Crabtree-positive and Crabtree-negative yeasts. We harmonized metabolites and reactions into a non-redundant shared ModelSEED ontological space, standardized compartmentalization, and refined gene-protein-reaction (GPR) rules. Using pathway-level visualization and systematic gap detection, we further improved the integrated network through literature-guided curation to correct stoichiometry, stereospecificity, and pathway architecture. Orthologous protein family reconstruction and functional annotation workflows were used to validate and inform GPR associations, with particular emphasis on ambiguous enzyme superfamilies and membrane-associated components. Using the resulting CFCMM, we built high-quality central carbon core models for each fungus and performed flux balance analysis to quantify ATP-yield variation under aerobic and anaerobic conditions, explicitly evaluating scenarios driven by differences in electron transport chain (ETC) composition. Simulations reproduced the expected fermentative yield of approximately 2 mmol ATP per mmol glucose under anaerobic conditions and separated the thirteen fungi into two bioenergetic groups under aerobic respiration based on Complex I status, with predicted yields of approximately 30 versus 22 mmol ATP per mmol glucose. Forcing flux through the alternative oxidase bypass further reduced ATP yields to approximately 12 and 4 mmol ATP per mmol glucose in Complex I-containing and Complex I-lacking fungi, respectively. Collectively, this work provides a manually curated, ModelSEED-consistent, and extensible fungal core metabolic template, deployed in DOE KBase as a resource for automated reconstruction of central carbon core models from any sequenced fungal genome. In addition, the CFCMM provides modular components for developing GEMs with more accurate energy predictions and enables robust comparative analyses of fungal bioenergetics and core metabolic diversity.

BackgroundThe epidermal growth factor receptor (EGFR) orchestrates highly context-dependent intracellular signaling networks whose architecture varies across cell types and is frequently rewired by targeted therapeutics. Systems-level reconstruction of these networks from phosphoproteomic data remains challenging because phosphorylation measurements identify signaling nodes but do not reveal the interaction paths that propagate signals between proteins. ResultsWe developed a computational framework integrating time-resolved phosphoproteomics with graph traversal algorithms to reconstruct EGFR-initiated signaling pathways across three contexts/conditions. A sign-assignment preprocessing procedure converts quantitative phosphorylation measurements into binary activation states across time points, defining a condition-specific active node set that filters the protein-protein interaction network. Breadth-First Search combined with interaction-weighted Beam Search is then applied to the STRING interaction database (v11.5) to enumerate candidate signaling paths. Applying this framework to phosphoproteomic datasets from EGF-stimulated HeLa cells, EGF-stimulated MDA-MB-468 triple-negative breast cancer (TNBC) cells, and EGF-stimulated MDA-MB-468 cells pretreated with the SHP2 inhibitor SHP099 yielded 260 paths in HeLa cells (117 unique topologies), 293 paths in MDA-MB-468 cells (155 unique), and 292 paths under SHP2 inhibition (85 unique). HeLa cells displayed a SRC-centered architecture dominated by ERBB2 and SHC1 first-hop effectors, converging on focal adhesion, HSP90 chaperone, CRKL adaptor, and integrin signaling arms. In contrast, MDA-MB-468 cells showed a PIK3CA/PTPN11 dual-axis architecture integrating direct PI3K engagement with SHP2-mediated GRB2-IRS1-ABL1 signaling. SHP2 inhibition abolished PTPN11-mediated pathways and induced PIK3CA dominance (69.2% first-hop), accompanied by compensatory ERBB3 engagement and a computationally predicted SYK/VAV1/LCP2 node set whose biological role warrants experimental validation. ConclusionsTime-resolved phosphoproteomics-guided BFS Beam Search over STRING interaction networks captures cell-type-specific EGFR signaling architectures and drug-induced pathway rewiring. This framework provides a systematic approach for transforming phosphoproteomic measurements into mechanistically interpretable signaling hypotheses specific to the cell-type-specific contexts, directly applicable to drug resistance modeling.

Difficult-to-treat (D2T) rheumatic disease affects approximately 12% of rheumatoid arthritis patients and resists sequential biologic therapy, yet no mechanistic model explains this refractoriness as a system-level phenomenon. Here we present the 3-Axis Integrative Framework (3-AIF), a six-variable ordinary differential equation system integrating mucosal tolerance, energy-gated neuroimmune danger sensing, and integrated stress response pathways coupled through Hill-function metabolic gating. Stability analysis reveals bistable dynamics with two co-existing attractors separated by a saddle point. Bifurcation analysis demonstrates fold catastrophe with hysteresis: recovery requires greater therapeutic effort than disease prevention. Sensitivity analysis identifies three dominant parameters mapping to neuroimmune activation, energy drain, and recovery capacity. Cross-disease transcriptomic consistency analysis across six public datasets (n=310, five rheumatic diseases, four tissue types) reveals compartment-specific axis dysregulation -- circulating cells show integrated stress response activation while target tissues show pathway exhaustion -- and disease-specific axis dominance patterns consistent with model predictions.

Spatial computational models such as agent-based models (ABMs) offer powerful in silico tools to study tumor dynamics, yet imaging data are still rarely used to inform these models directly. We present an ABM optimization framework that leverages convolutional encoders to compare spatial patterns between experimental imaging data and ABM-generated outputs within a shared latent space. This quantitative comparison was used to estimate ABM parameters across three datasets, ranging from synthetic data to 3D tumoroid-T cell co-culture microscopy and histopathology images from The Cancer Genome Atlas skin cutaneous melanoma samples. Estimated parameters were evaluated using data-derived features and experimental knowledge, including experimental conditions and gene expressions. Simulations using optimized parameters reproduced key spatial features of the training images, such as tumor boundary complexity and tumor-tumor neighborhood structure. Together, these results demonstrate a flexible framework for ABM parameter optimization using spatial data across modalities, enabling systematic investigation of how spatial architecture influences tumor progression and immune interactions.

Protein tyrosine phosphatase receptor J (PTPRJ) restrains cell proliferation and migration by dephosphorylating receptor tyrosine kinases (RTKs) including the epidermal growth factor receptor (EGFR). PTPRJ is a purported tumor suppressor, and alterations to its expression and/or function are associated with colorectal, breast, lung, and other cancers. While there is interest in controlling PTPRJ-regulated phenotypes, efforts are limited by the complexity of PTPRJ-mediated signaling. PTPRJ dephosphorylates multiple RTKs, and the degree to which PTPRJ control of signaling and phenotypes depends on local cellular RTK activation profiles is unknown. To probe the context dependence of PTPRJ signaling regulation, we collected signaling measurements across 16 pathway nodes at two time points in a panel of HSC3 carcinoma cells engineered with different PTPRJ expression profiles. Cells were treated with three different RTK ligands, and paired phenotype measurements (viability, wound healing, xCELLigence cell index) were made. Partial least squares regression models were developed to predict relationships between PTPRJ-regulated signaling pathways and cell phenotypes. The model effectively separated contributions to variance arising from the PTPRJ expression background and growth factor context. In testing model predictions, we demonstrated that PTPRJ suppressed MET-induced cell cell proliferation via regulation of a HER3/AKT signaling axis that stabilized PTPRJ expression through an unanticipated feedback mechanism. We also found that PTPRJ regulated HSC3 cell migration via JNK signaling that was preferentially activated by MET. Our results identify new regulatory nodes through which PTPRJ influences cancer cell phenotypes and demonstrates that these processes preferentially occur in the context of distinct RTK activation states.

BackgroundThe complexity of gene regulatory networks, involving thousands of genes, poses a fundamental challenge to understanding cancer phenotype reversal. However, recent evidence suggests that the effective dimensionality of normal and tumor transcriptional manifolds is low, and that small panels of genes can discriminate perfectly between normal and tumor samples. MethodsWe build upon two previously developed concepts: (i) highly accurate normal and tumor gene markers (namely, N-and T-markers), defined as genes with exclusive expression intervals in normal and tumor samples, respectively; and (ii) gene deregulation networks (GDNs), represented as directed acyclic graphs encoding causal relationships between gene deregulation events. A subset of genes appearing in both marker classes (NT-markers) act as bridging nodes between the N-and T-GDNs. Starting from these elements, we introduce a quantitative dynamical model based on node frequency and connectivity to assess how gene intervention effects propagate through the GDN and thereby predict their overall impact on the tumor tissue. ResultsAccording to the model, interventions on pure T-markers (T-markers that are not NT-markers) produce effects largely confined to the T-GDN, with a minimal perturbation of the N-network. Interventions on pure N-markers (N-markers that are not NT-markers) generate a perturbation of both networks, but with limited effect. In contrast, interventions on NT-markers with high activation frequency in both tumors and normal state (e.g., AGER in lung adenocarcinoma: 98% in tumor samples, 75% in normal samples) can induce bidirectional phenotype shifts. For an effective combination of targets, coverage across tumor samples must be maximized. At the same time, in the T-GDN the number of nodes unreached by the reverse cascade following the intervention must be minimized, as these regions may act as escape routes for the tumor. Escape probability further depends on the tumor stage and the tumors activation rate of new T-genes. When targeting NT genes, high frequency in normal samples should also be prioritized. ConclusionsHigh-frequency NT-genes, due to dual network connectivity and tissue relevance, represent optimal targets for achieving at least partial phenotype reversal. This framework provides a quantitative guide for prioritizing gene therapy targets and designing combination strategies that balance coverage, escape minimization, and normal tissue relevance.

Agent-based modeling (ABM) is a computational method for predicting the emergent outcomes of interacting, autonomous individuals in a complex system. Here, ABM is used to simulate interactions between fibroblast and myofibroblast cells during idiopathic pulmonary fibrosis (IPF) in alveolar tissue microenvironments. These microenvironments are derived from histology of a healthy human lung sample and moderate- and severe-IPF lung samples. Fibroblast differentiation, cell migration, and collagen secretion in response to the spatial distribution of the cytokine transforming growth factor-beta are captured in the ABM using NetLogo software. Results are presented from one simulated year without treatment and with mechanisms representing treatment by pirfenidone and pentoxifylline, alone and in combination. A total of 180 in silico experiments are run, analyzed, and compared in a high-throughput workflow. The effects of the initial number of fibroblasts and treatment scenarios on various metrics related to collagen accumulation and collagen invasion into alveolar regions are determined. The ABM and the analysis files are shared to facilitate model reuse. By integrating computational modeling of IPF and therapeutics, this research aims to improve understanding of fibrosis progression and assess the efficacy of novel and existing treatments targeting different mechanisms to inform decision-making for IPF treatment.

Stachys affinis (Chinese artichoke) tubers contain 50-80% stachyose by dry weight, the most concentrated dietary source of raffinose-family oligosaccharides (RFOs) known. Because humans lack sufficient -galactosidase activity, stachyose transits intact to the colon where microbial fermentation yields short-chain fatty acids (SCFAs). However, the quantitative impact of stachyose-derived SCFAs on host hepatic energy metabolism has not been systematically explored using genome-scale metabolic models. Three stachyose dose scenarios (Low/Mid/High: [~]25, 50, 100 g fresh tubers) were translated to SCFA availability vectors. Hepatic metabolic responses were simulated using Recon3D (10,600 reactions) and Human-GEM (13,417 reactions) under strict hepatocyte-like media, maximizing ATP maintenance flux (ATPM). FVA across multiple optimality thresholds (90-100%) and pFBA confirmed solution robustness. One-at-a-time sensitivity analysis characterized ATPM responses to individual parameter perturbations, and a ratio sensitivity sweep across six alternative SCFA profiles assessed dependence on assumed fermentation ratios. A targeted rescue experiment addressed model-specific propionate catabolism gaps. Both models showed dose-dependent ATPM increases (Recon3D: +71 to +286%; Human-GEM: +103 to +413% above baseline), with the 19-33% inter-model gap attributable entirely to Human-GEMs functional propionate catabolism pathway. FVA confirmed near-unique optimal solutions (ATPM ranges [~]1% at 99% optimality, widening to [~]10% at 90%). Parsimonious FBA preserved identical ATPM values while reducing total flux by [~]4-14%, confirming objective robustness. SCFA ratio sensitivity across six alternative profiles showed 27- 28% ATPM variation, indicating qualitative robustness. Butyrate yielded the highest ATP per mole ([~]22) in both models; propionate sensitivity was zero in Recon3D but [~]15.25 mmol ATPM/mmol propionate in Human-GEM. Reopening propionyl-CoA carboxylase (PPCOACm) in Recon3D under strict constraints converged ATPM to within 0.3-0.7% of Human-GEM, cross-validating both reconstructions. This reproducible dual-model pipeline identifies model-specific pathway gaps and provides cross-validated predictions to guide future experimental studies of how dietary SCFAs influence hepatic ATP metabolism.

SummaryBoolean and logic-based modeling approaches are well suited for the analysis of complex biological systems, particularly when detailed biochemical and kinetic information is unavailable. In such settings, biological pathways are represented as networks capturing system components and their interactions, providing a simplified yet informative abstraction of system behavior. While the structural topology of these networks is often well characterized, the absence of mechanistic detail limits the applicability of parameter-dependent modeling frameworks. To address this, we present BoolDog, a Python package for the construction, simulation, and analysis of Boolean and semi-quantitative Boolean networks. BoolDog supports synchronous simulation with events, attractor and steady-state identification, network visualization, and the systematic transformation of logic-based models into continuous ordinary differential equation (ODE) systems -- enabling the seamless integration of discrete and continuous modeling paradigms. Networks can be imported and exported across standard formats, and BoolDog integrates natively with established Python libraries for network analysis and visualisation, including NetworkX, igraph, and py4Cytoscape. Together, these capabilities provide a flexible, accessible, and interoperable platform for logic-based modeling of complex biological systems. Availability and implementationBoolDog is implemented in Python and available at https://github.com/NIB-SI/BoolDog/.

Quantitative Systems Pharmacology (QSP) models play an inherently interventional role in pharmaceutical research and development, functioning as executable causal systems for designing, evaluating, and replacing clinical trials. However, deploying QSP as an experimental planning engine remains constrained by the difficulty of translating unstructured literature descriptions of clinical or preclinical scenarios into reproducible, simulation-ready model interventions. Motivated by this issue, we propose an agent-based framework that operationalizes QSP models as intervention-ready experimental systems by automatically extracting and executing literature-derived scenarios. The framework combines semantic grounding of model entities with a large language model (LLM)-driven Scenario Extractor and a dual-agent Scenario Mapper. Rather than relying on opaque, single-shot reasoning, our pipeline converts free-text interventions into precise parameter configurations through discrete, verifiable work orders. Moreover, our dynamic Human-in-the-Loop (HITL) strategy empowers modelers to resolve biological ambiguities interactively. Across four diverse kinetic ordinary differential equation (ODE)/QSP models and seven Subject Matter Expert (SME)-curated literature scenarios, our model resolved all selected scenarios into correct executable parameter changes, including multi-dose interventions, unit conversions, no-op scenarios, and ambiguity-triggered HITL cases, demonstrating that structured collaboration between experts and agentic systems can resolve scenarios that standalone raw Systems Biology Markup Language (SBML) reasoning LLM calls handle unreliably.